##########################################################################################

library(data.table)
library(optparse)
library(Seurat)
library(monocle)
library(ArchR)
library(GenomicFeatures)
library(dplyr)

##########################################################################################
option_list <- list(
    make_option(c("--rna_data_file"), type = "character"),
    make_option(c("--rna_data_monocle_file"), type = "character"),
    make_option(c("--atac_data_file"), type = "character"),
    make_option(c("--time_diff_gene_file"), type = "character"),
    make_option(c("--gff3_file"), type = "character"),
    make_option(c("--gene_ensg_file"), type = "character"),
    make_option(c("--out_path"), type = "character") 
)

if(1!=1){
    
    rna_data_file <- "~/20231121_singleMuti/input/testis_combined.annotationCellType.Rdata"

    ## 单细胞表达文件
    rna_data_monocle_file <- "~/20231121_singleMuti/results/monocole/testis.monocle.Rdata"

    ## atac表达文件
    atac_data_file <- "~/20231121_singleMuti/results/atac_res/testis_combined_peak.Rdata"

    ## 差异表达文件
    time_diff_gene_file <- "~/20231121_singleMuti/results/monocole/pseudotime_ordergene.tsv"

    ## 注释文件
    gff3_file <- "~/ref/GTF/gencode.v32.annotation.gff3"
    gene_ensg_file <- "~/ref/GTF/gencode.v32.gene_ensg.tsv"

    ## 输出
    out_path <- "~/20231121_singleMuti/results/monocole"

}

###########################################################################################
parseobj <- OptionParser(option_list=option_list, usage = "usage: Rscript %prog [options]")
opt <- parse_args(parseobj)
print(opt)

rna_data_file <- opt$rna_data_file
rna_data_monocle_file <- opt$rna_data_monocle_file
atac_data_file <- opt$atac_data_file
time_diff_gene_file <- opt$time_diff_gene_file
gff3_file <- opt$gff3_file
gene_ensg_file <- opt$gene_ensg_file
out_path <- opt$out_path

dir.create(out_path , recursive = T)

###########################################################################################

a <- load(rna_data_file)
## scrnat
b <- load(atac_data_file)
## testis_combined_peak
c <- load(rna_data_monocle_file)
## ALL_cds_all

## 注释文件
txdb <- makeTxDbFromGFF(gff3_file)
gene_ensg <- read.table(gene_ensg_file)

## 逆时序差异基因列表
Time_diff <- data.frame(fread(time_diff_gene_file))

###########################################################################################
## RNA的时间
cell_time <- data.frame(cell_id = ALL_cds_all$cell , Pseudotime = ALL_cds_all$Pseudotime )
cell_time$cell_id <- gsub( '_' , '#' , cell_time$cell_id)
## 时间放到100之内,atac需要
cell_time$Pseudotime <- cell_time$Pseudotime/(max(cell_time$Pseudotime)) * 100

################################################################################
## 构建基因的位置信息
gene_gr <- genes(txdb, columns = c("GENEID"))
colnames(gene_ensg) <- c("GENEID" , "name")

gr_df <- data.frame(
  seqnames = seqnames(gene_gr),
  start = start(gene_gr),
  end = end(gene_gr),
  strand = strand(gene_gr),
  GENEID = names(mcols(gene_gr)$GENEID)
)

## gene name对应多个ensg的，选第一个
gr_df <- merge( gr_df , gene_ensg , by = "GENEID" )
gr_df <- gr_df %>% 
group_by( name ) %>%
summarize( seqnames = seqnames[1] , start = start[1] , end = end[1] , strand = "*"  )

gr <- GRanges(
  seqnames = Rle(gr_df$seqnames),
  ranges = IRanges(start = gr_df$start, end = gr_df$end),
  strand = gr_df$strand,
  name = gr_df$name
)

names(gr) <- gr_df$name

################################################################################

use_gene <- rownames(scrnat@assays$RNA$counts)[rownames(scrnat@assays$RNA$counts) %in% names(gr)]
gr <- gr[use_gene]

seruat_cds_all <- SummarizedExperiment(scrnat@assays$RNA$counts[use_gene,] , rowRanges = gr)
colnames(seruat_cds_all) <- gsub( "_" , "#" ,  colnames(seruat_cds_all) )

## 整合atac和表达，按照每个细胞合并
proj <- addGeneExpressionMatrix(input = testis_combined_peak, seRNA = seruat_cds_all, force = TRUE)
projHeme3 <- addImputeWeights(proj)

################################################################################
## 拟时序，rna的时序
projHeme5 <- projHeme3
projHeme5@cellColData[cell_time$cell_id,"Pseudotime"] <- cell_time$Pseudotime

################################################################################

out_file <- paste0( out_path , "/testis_combined_peak.combineRNA.Rdata" ) 
testis_combined_peak_combineRNA <- projHeme3
save( testis_combined_peak_combineRNA, file = out_file )

